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ABSTRACT 


Invasion, spread and establishment of invasive alien species in a new environment 
causes serious ecosystem perturbation and native species extinctions. The problem is 
further aggravated under climate change as some invasive alien species perform better 
under elevated temperature and carbon dioxide regimes. Currently unsuitable regions 
such as high-altitude areas and mountains are likely to become more suited to invasion 
under future climatic conditions. We have modelled the distribution of one of the most 
invasive alien species, parthenium weed (Parthenium hysterophorus L.) which is rapidly 
colonizing different parts of Bhutan. The study was implemented in R environment using 
BIOMOD2 package; an ensemble modelling platform for species distribution modelling. 
Under current climate scenario, about 2.83% (1,099.01 km?) of the country’s total area 
was predicted to be suitable for parthenium weed invasion, covering 17 out of the 20 
districts in Bhutan. Under future climate scenarios, the highest suitability was predicted 
under RCP4.5 2050 period with about 5,419.69 km? anticipated to be suitable. Except for 
Bumthang, all districts showed suitability to invasions under future climate scenarios. 
Generally, districts located in the west and south showed more suitability than those in 
the east and central region. The highest elevation of suitability was predicted to be at 
2,931 m above sea level; an upward shift of about 753 m. Based on these findings, there 
is an urgent need to develop management programs and raise public awareness on the 
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adverse impacts of parthenium weed in Bhutan. 


Introduction 


The introduction of alien species outside of their na- 
tive range deliberately or accidentally destabilizes eco- 
system functioning [1]. Native species extinctions have 
been well-documented [2] and debate on whether alien 
species invasion could lead to mass extinctions of native 
species is growing [3]. Even the Intergovernmental Plat- 
form for Biodiversity and Ecosystem Services (IPBES) of 
the United Nations has recognized the threats posed by 
biotic invaders [4]. Further, climate change is escalating 
the situation by removing climate barriers [5], facilitat- 
ing the spread and establishment [6] and shifting ranges 
[7]. The second most threat after habitat loss to biodi- 
versity is Invasive Alien Species (IAS) [8]. 


Bhutan, a tiny Himalayan democratic constitutional 
monarchy with an area of 38,394 km’, a population of 
about 0.72 million [9] endowed with very rich biodiver- 
sity is no exception and considered spared from inva- 
sion by JAS. There is already a number of IAS reported in 
Bhutan [10-16], including Crofton weed (Ageretina ade- 
nophora (Spreng.) King & H.Rob.), Siam weeds (Chromo- 
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laena odorata L.) R.M.King & H.Rob.), lantana (Lantana 
camara L.), mile-a-minute (Mikania micrantha Kunth) 
and parthenium weed (Parthenium hysterophorus L.). 
Between 1974 to 2000, from a total of 127 species im- 
ported into the country, four species of trees and shrubs 
and 18 species of grasses have become very serious 
weeds [15]. There are already about 43 Invasive Alien 
Plant Species (IAPS) recorded in Bhutan [17] and 16 are 
a major IAPS [16]. The most comprehensive inventory 
of IAPS undertaken by [10] reported about 101 species 
of IAPS in the country. 


The detrimental effects of IAPS to public health, agri- 
culture, fisheries, and the economy is well-known with 
losses estimated in millions of dollars. For instance, IAPS 
such as parthenium weed is known to cause allergies, 
asthma and bronchitis, and in some cases even death 
[18]. Consumption of Crofton weed has led to the death 
of 30 horses [19] and farmers in Bhutan have reported 
impacts of IAPS on highland pastures [20] and on crops 
and forests [21]. In Australia, it is estimated that inva- 
sive weeds alone cause crop production loss to the tune 
of AU$ 1.271 billion per year [22], while in Nepal it is 
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estimated that IAPS cause losses about US$ 1.4 billion per 
year [23] and India, about US$ 91 billion per year [24]. 


One of the most invasive alien species known to the world 
with its origin in central Mexico [25] and reported in 45 
countries, including Bhutan [26,27] is the parthenium 
weed. In Bhutan, it was first reported in 1992 along road- 
sides, in crop fields and settlement areas [28] and it is 
now reported at elevation as high as 2,600 m Above Sea 
Level (ASL) [13]. To locals in eastern Bhutan, it is known 
as “Amartala ngon” (Weed from the place in India called 
Amartala/Agartala) or “Gungthri ngon’” (seeds fallen from 
the sky) (S. Wangdi, personal communication, 12 March 
2021). It reduces crop yield by as much as 50% [27] and 
causes detrimental effects to public health, agriculture, 
forestry, and livestock. When used as a livestock feed, milk 
becomes bitter [29] and meat can become tainted. Under 
climate change, it is likely to spread because of the ben- 
efit accrued from global warming and atmospheric CO2 
enrichment [30]. Many invasive species can expand their 
range under climate change [31,23]. Therefore, there is a 
need to accurately predict the potential distribution and 
understand the impact of parthenium weed in Bhutan. 
This will not only help to make informed decisions on 
matter concerning biodiversity, public health, agriculture 
and the economy but also help in early detection to enable 
prompt actions in order to reduce management costs, 
which otherwise would become expensive or impossible 
at later stages. 


Other than the study conducted by [32] and [33], no study 
had been conducted on IAPS distribution in Bhutan. Both 
the studies implemented the most commonly used spe- 
cies distribution modelling program (-MaxEnt) [34] and 
provided a general overview of parthenium weed distri- 
bution at the national level and West-Central region, re- 
spectively under current and a future climate (predicted 
for 2070). To develop management programs and inform 
policy-decision, further details such as areas under high 
risk of invasion or currently invaded, percent invasion 
under current climate and future climates, elevational 
shifts and distribution predicted for short- (2030s) and 
medium-term (2050) periods is necessary. Using the BIO- 
MOD2 package [35,36] and ensemble approach in R v 
4.0.3 [37], we have modelled the distribution of partheni- 
um weed in Bhutan for the current (1979-2013) and two 
future climate periods (2050 and 2070). Additionally, we 
have analyzed the distribution of parthenium weed at the 
district level, which will be important while planning for 
early detection, prevention and control and management 
programs. 


Materials and Methods 


Spatial thinning of occurrence points 


The occurrence points of parthenium weed were collect- 
ed from the three national highways: Thimphu to Sarpang 


national highway, Thimphu to Punakha national highway 
and Thimphu to Samtse national highway via the Phuent- 
sholing to Samtse internal road (Figure 1). The occur- 
rence points were collected using a handheld GPS with 
3 m accuracy. The total distance surveyed was about 360 
km. The records were also taken for the dried and dead 
remains of the parthenium weed from the previous years 
in case the data collection time coincided with the early 
emergence time at that site. A total of 151 occurrence 
points were collected between April to August 2020. As 
spatial autocorrelation can lead to model over fitting and 
can reduce model performance, we tested the data for 
spatial autocorrelation and thinned in ArcGIS 10.4 using 
SDMTOOLBOX v2.4 [38]. A total of 55 occurrence points, 
which were retained after spatial thinning from 151 re- 
cords were finally used for modelling. 
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Figure 1. Occurrence points (red dots) along the surveyed 
roads (green). Note: ( « ) Occurrence, (—) surveyed roads, ( 
—) Highway (~~) Districts 


Selection of predictor variables and testing for 
multicollinearity 


The likely distribution of a species is commonly predict- 
ed by using bioclimatic variables and the WORLDCLIM 
2 (http://worldclim.org) [39] is the source mostly used 
to obtain these variables. We have used the CHELSA (Cli- 
matologies at High Resolution for the Earth’s Land Sur- 
face Areas) website (https://chelsa-climate.org/) [40] to 
obtain the 19 bioclimatic variables; only being that it has 
out-performed the WORLDCLIM 2 in terms of precipita- 
tion, which was demonstrated from a preliminary study 
using Bhutan [40]. The 19 bioclimatic variables at 30 arc 
sec (ca. 1 km’) resolution were downloaded for the cur- 
rent (1979-2013), and two future climate periods, e.g. for 
2050 (2041-2060 average) and 2070 (2061-2080 aver- 
age) for Representative Concentration Pathway (RCP) 4.5 
and RCP8.5. 


The RCP4.5 assumes a steady increase in radiative forcing 
with a projected global mean surface temperature of 1.4°C 
to 1.8°C, while the RCP8.5 assumes the highest increase in 
radiative forcing with projected global mean surface tem- 
perature of 2.0°C to 3.7°C [41]. The choice of RCPs was 
also in line with the climate change projection for Bhutan 
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[42]. Given that Bhutan is a mountainous country with 
high microclimate variability, we have used the global cir- 
culation model, MIROC5 has used by others working in 
the Himalayan region [43-47]. Additionally, we have gen- 
erated two topographic variables (slope and aspect) from 
the Digital Elevation Model (DEM), which was obtained 
from the Shuttle Radar Topographic Mission (SRTM). 


Predictor variables were tested for multicollinearity [48] 
using the vifcor function in the usdm package [49] exe- 
cuted in R v 4.0.3 [37]. There is no consensus as to what 
an acceptable threshold is, and how it could be used for 
discarding correlated variables. Based on the recommen- 
dation by [50] we have used a threshold of r (Pearson’s 
correlation) >0.8 to screen out the correlated variables. 
The variable retained after multicollinearity test were: 
(i) Mean Diurnal Range (Bhutan_Bio2), (ii) Isothermality 
(Bhutan_Bio3), (iii) Min Temperature of Coldest Month 
(Bhutan_Bio6), (iv) Mean Temperature of Wettest Quar- 
ter (Bhutan_Bio8), (v) Precipitation Seasonality (Bhu- 
tan_Bio15), (vi) Precipitation of Warmest Quarter (Bhu- 
tan_Bio18), (vii) Precipitation of the Coldest Quarter 
(Bhutan_Bio19), (viii) Aspect, and (ix) Slope. The study 
did not include other anthropogenic variables such as 
roads, land use and settlements as future projections of 
these variables were not available. 


Modelling approach and calibration 


There are many species distribution models used to 
model species distributions in a geographic space and 
environment [51]. We have used an ensemble model- 
ling approach, as it accounts for uncertainties of differ- 
ent models and reported to perform better than a single 
model [52,36]. The package BIOMOD (BIOMOD2 package 
in R) developed by [35] has been especially designed to 
address this issue. The package provides 10 models but 
for our study, we have only used five, viz., General Additive 
Model (GAM), Generalized Linear Model (GLM), General- 
ized Boosted Models (GBM), maximum entropy (MAX- 
ENT:Phillips) and Random Forest (RF). These models are 
the most used and are known to perform well [53,54]. 


The above models in BIOMOD2 were run with these set- 
tings. While formatting, the number of pseudo-absences 
generation was repeated three times, with their number 
generated set to 2000. Keeping a buffer of 10 km from the 
presence points, pseudo-absences were randomly select- 
ed [55]. Pseudo-absences were required as the models 
used were based on presence-absence data and we did 
not have absence data. At the modelling step, we calibrat- 
ed the models by splitting the data 70:30 (70% for model 
calibration/training and 30% for model evaluation /test- 
ing) [56]. In total, we built 75 models (three pseudo-ab- 
sence runs, five models and five evaluation runs). In order 
to build ensemble models, models below True Skills Sta- 
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tistics (TSS) value of 0.6 and Area Under the Curve (AUC) 
0.8 were screened out and weights were assigned to each 
model balanced to their evaluation metric score; what’s 
called a weighted-mean approach. Using TSS scores as a 
cut-off value, we produced binary maps of suitable and 
unsuitable areas. Change in suitability areas under future 
climatic conditions was produced by using BIOMOD2 
range size function. All the maps produced in R environ- 
ment using BIOMOD2 were imported into ArcGIS 10.4 for 
analysis and visualization. 


Model evaluation 


We have used the AUC of the Receiver Operating Charac- 
teristics (ROC) [57] and TSS to evaluate the model per- 
formance. The AUC is a measure of the ratio between the 
observed presence-absence values and the model predic- 
tions and the values range from 0 to 1. This is a threshold 
independent measure and according to the scores, a mod- 
el is classified as excellent (0.9-1.0), good (0.8-0.9), fair 
(0.7-0.8), and poor (0.6-0.7) [58,59]. TSS is a threshold 
dependent measure and its values range from -1 to +1. 
Based on TSS, a model is categorized as no better than 
random (< = 0), poor (0.2-0.5), useful (0.6-0.8), excellent 
(>0.8) and prefect agreement (+1) [60,61]. 


Results 


Prediction accuracy of models and variable impor- 
tance 


The score of individual models based on average AUC 
ranged from 0.86 (MAXENT:Phillips) to 0.97 (GBM/RF); 
while the average TSS scores ranged from 0.73 (MAXENT. 
Phillips) to 0.92 (GBM) (Table 1). Based on the AUC score, 
the best performing model was the GBM and the RF while 
based on the TSS score, it was GBM. Nonetheless, all the 
single models performed better than random models but 
not better than an ensemble model (Figure 2). The medi- 
an AUC of the ensemble model was 0.99 while TSS was 
0.97. As predicted by the model, compared to the current 
period, the suitable habitat area of parthenium weed 
will increase by the years 2050 and 2070. Amongst the 
nine predictor variables used for modelling, mean tem- 
perature of the wettest quarter (Bhutan_Bio8) was the 
most important variable influencing parthenium weed 
distribution in Bhutan. Precipitation seasonality (Bhu- 
tan_Bio15), precipitation of the warmest quarter (Bhu- 
tan_Bio18) and the min temperature of the coldest month 
(Bhutan_Bio6) influenced this to some extent (Table 2). 
Parthenium weed suitability increased with increasing 
mean temperature of wettest quarter (Bhutan_Bio8) and 
precipitation seasonality (Bhutan_Bio15) while it de- 
creased with increasing precipitation of warmest quarter 
(Bhutan_Bio18) and the min temperature of the coldest 
month (Bhutan_Bio6) (Figure 3). 
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Table 1. Performance score of individual model based on the Receiver Operating Characteristics (ROC) and the Area Under 
the Curve (AUC) of and the True Skills Statistics (TSS). AUC and TSS presented is the average of five evaluation runs and three 
pseudo-absences run of each modelling algorithm. 


Model AUC TSS 
GBM 0.97 0.92 
RF 0.97 0.89 
GLM 0.95 0.88 
GAM 0.92 0.82 
MAXENTPhillips 0.86 0.73 
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Figure 2. Performance of individual model and ensemble model based on the Receiver Operating Characteristics (ROC) and 
the Area Under Curve (AUC) and True Skills Statistics (TSS) in predicting the distribution of P. hysterophorus in Bhutan. 
Models used: Generalized Additive Model (GAM), Generalized Linear Model (GLM), Generalized Boosted Model (GBM), Ran- 
dom Forest (RF), and Maximum Entropy (MAXENT-Phillips). The bottom of the whisker of the box-plot represents minimum, 
bottom edge represents lower quartile, bold line inside represents median, upper edge represents upper quartile, top of the 
whisker represents maximum and filled dots represents outliers. The predictive performance of the ensemble model (AUC = 
0.99, TSS = 0.97) is represented by the coloured dashed line. 


Table 2. Importance of predictor variables used in distribution modelling of P. hysterophorus in Bhutan. Variable importance 
is presented as the mean over five evaluation runs and three pseudo-absences datasets of each modelling algorithm, and as 
the mean of means amongst them. The models: Generalized Additive Model (GAM), Generalized Boosted Model (GBM), Gen- 
eralized Linear Model (GLM), Random Forest (RF), Maximum Entropy (MAXENT-Phillips). Variables of greatest influence are 
highlighted in bold. 


Model Aspect * Bio15 | *_Bio18 | *_Bio19 | *_Bio2 * Bio3 * Bio6 * Bio8 Slope 


GAM 0.1353 0.5723 0.7856 0.2431 0.2253 0.2691 0.3478 0.8304 0.0675 
GBM 0.0285 0.0993 0.0142 0.2921 0.0320 0.0039 0.0069 0.6434 0.0027 
GLM 0.1277 0.4938 0.3348 0.1929 0.1383 0.1040 0.3701 0.9205 0.0199 
RF 0.0149 0.0379 0.0240 0.1437 0.0352 0.0097 0.0701 0.2038 0.0039 


MAXENT. | 0.1937 | 0.2331 0.4140 |0.1498 | 0.1659 0.0543 0.6685 0.0685 0.1455 
Philipps 


Mean of | 0.1000 0.2873 0.3145 0.2043 0.1193 0.0882 0.2927 0.5333 0.0479 
means 


Note: *Bhutan 
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Figure 3. Response curve of probability of occurrence obtained for top four contributing predictor variables of ensemble 
model by weighted-mean approach based on Area Under Curve (AUC). Precipitation seasonality (Bhutan_Bio15), precipita- 


tion of warmest quarter (Bhutan_Bio18), min temperature of 
wettest quarter (Bhutan_Bio8). 


Current and future distribution of parthenium 
weed in Bhutan 


The potential distribution of parthenium weed in Bhutan 
based on the current climatic conditions and occurrences 
points showed that about 1,099.1 km? (2.83%) of the total 
area of Bhutan is suitable for parthenium weed invasion 
(Figure 4 and Table 3). This includes regions within 17 
out of the 20 districts of Bhutan. Amongst districts, Wang- 
duephodrang had the highest suitability of about 302.09 
km? (27.49%) followed by Punakha, Dagana and Samtse 
with 142.68 km? (12.98%), 129.31 km? (11.76%) and 
107.01 km? (9.74%), respectively. Gasa and Haa had the 
lowest suitability. 


Under future climatic projections, the land suitable for 
parthenium weed invasion increased in both the climate 
periods 2050 and 2070 (Figure 4 and Table 3). The high- 
est suitability was predicted under RCP4.5 2050 with 
about 5,419.69 km? (13.95% of Bhutan’s area), followed 
by RCP8.5 2070 with about 3,954.97 km? (13.01% of 
Bhutan’s area). The lowest is predicted under RCP4.5 
2070 with about 1,771.27 km? (Table 3). There was also 
a change in the number of districts suitable to parthe- 
nium weed invasion. The 17 districts that were suitable 
under current remained stable under RCP4.5 2050 sce- 
nario while additional two districts (Paro and Thimphu) 
became suitable under RCP8.5 2050. However, under 
RCP4.5 2070 scenario, only one additional district (Paro) 
became suitable compared to current prediction while 
under RCP8.5 2070, districts’ suitability remained same 
as under RCP8.5 2050 scenario. Bumthang remained free 
from parthenium weed invasion under any period. 


District-wise, under RCP4.5 2050 scenario, Sarpang had 
the highest suitability with about 706.72 km? followed 
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the coldest month (Bhutan_Bio6), and mean temperature of 


by Dagana with about 692.23 km?. Under RCP8.5 2050 
scenario, Dagana had the highest suitability with about 
462.60 km? followed by Wangduephodrang with about 
429 km*. Under RCP4.5 2070 scenario, Dagana had the 
highest suitability with about 441.42 km? followed by 
Wangduephodrang with about 434.74 km?. Under RCP8.5 
2070 scenario, the highest suitability was predicted in Da- 
gana with about 681.08 km? followed by Wangduepho- 
drang with about 560.69 km’. 


In terms of change in suitability between current and fu- 
ture climate periods, the highest gain in suitability was 
determined under RCP4.5 2050 scenario with about 
4,391.94 km? while the lowest was determined under 
RCP 4.5 2070 scenario with about 1,959.65 km? (Figure 4 
and Table 4). In terms of loss, the highest loss in suitability 
was determined under RCP4.5 2070 with about 178.35 
km? while the least was determined under RCP4.5 2050 
(Figure 4 and Table 4). Similarly, the highest stability was 
predicted under RCP4.5 2050 with about 1,130.31 km? 
while the lowest was predicted under RCP4.5 2070 with 
about 957.53 km? (Figure 4 and Table 4). 


The difference in suitability in terms of elevation range at 
different climate periods was also determined. There was 
a stark difference in the elevation range of parthenium 
weed distribution between current and future climates. 
The occurrence points were collected within the elevation 
range between 210-1,423 m asl but ensemble model pre- 
dicted the elevation range for current distribution to be 
between 109-2,178 m asl. The upper elevation range in- 
creased under both the future climates and scenarios and 
ranged between 2,513-2,931 masl. The highest elevation 
of about 2,931 m asl was predicted under RCP8.5 2070; 
an upward migration by about 753 m asl as compared to 
the current distribution. 
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Figure 4. Ensemble maps showing current and future distribution of P. hysterophorus in Bhutan. Leaf green areas indicate 
suitable areas while white areas indicate not suitable areas. The maps were prepared in ArcGIS 10.4. Note: (@™@) Suitable, ( 


) Not suitable, (4) Districts 


Table 3. Current and future predicted habitat suitability areas in Bhutan. 


Period Suitability | Area change | Percent (%) of current | Percent (%) of Status 
(Km?) (Km?) predicted distribution | Bhutan’s area 

Current 1099.10 2.83 

RCP 4.5 2050 5419.70 4320.60 393.10 13.95 Gain 

RCP 8.5 2050 3364.18 2265.08 206.09 8.66 Gain 

RCP 4.5 2070 2870.37 1771.27 161.16 7.39 Gain 

RCP 8.5 2070 5054.07 3954.97 359.84 13.01 Gain 


Table 4. Change in the distribution of P. hysterophorus between current and future periods. 


Period Loss Stable Gain 
RCP 4.5 2050 5.57 1130.31 4391.94 
RCP 8.5 2050 27.87 1108.02 2334.19 
RCP 4.5 2070 178.35 957.53 1959.65 
RCP 8.5 2070 82.49 1053.40 4051.95 
Discussion [63,64], number of occurrences used for model training 


While species distribution modelling is a popular choice 
to assess potential risk of species invasion [62], there 
are various factors that influence the performance and 
thereby the outcome of the models. These include the ex- 
tent and resolution of the study area and threshold used 


[65], species type - whether common or rare [66], vari- 
able selection technique, environmental predictors and 
species distributional data [67], choice of models [68], 
modelling approaches [69], GCMs used [70], model com- 
plexity [71], collinearity [72], presence-only or presence/ 
absence data used [62,73], method and number of back- 
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ground points used [55,62], method of evaluation [74] 
and binary maps’ thresholds used [75]. Several algorithms 
are recommended for modelling as no single model per- 
forms best under all situations [76]. Also, ensemble mod- 
els account for uncertainties in predictions of different 
algorithms and performs better than a single modelling 
algorithm [52,36]. Therefore, we decided to use ensemble 
modelling approach in our study. We used both the AUC 
and TSS to evaluate models as the use of single AUC to 
pick a model is not recommended [60]. The performance 
of individual model ranged from good to excellent based 
on the average AUC, which ranged from 0.86 to 0.97. True 
skills statistics ranged from 0.73 to 0.92 indicating useful 
to excellent performance. The prediction accuracy was 
significantly improved by combining models and building 
an ensemble in which the AUC was 0.99 (excellent) and 
TSS was 0.97 (excellent). 


Given that Bhutan values its biodiversity conservation 
and environmental protection, ensuring both are main- 
streamed into the country’s development programs and 
policies, the use of species distribution modelling seems 
inevitable. Yet, species distribution modelling has nev- 
er been used, especially in spatial planning and deci- 
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sion-making in the context of plant invasion [12] tried to 
model six invasive alien species, including parthenium 
weed in Bhutan using the most commonly implement- 
ed MaxEnt tool [34]. However, the study did not provide 
appropriate details in terms of districts suitability and 
elevational changes other than range expansion, contrac- 
tion, and directional change, which have limited use when 
planning for management programs. Similarly, using the 
MaxEnt tool, [33] studied parthenium weed distribution 
and its impact in the West-Central region but did not 
provide much detail. Also, they considered only one fu- 
ture period (2070). We have modelled for the first time 
the distribution of parthenium weed in Bhutan based on 
ensemble modelling approach using the BIOMOD2 pack- 
age developed by [35] and [36] in R v 4.0.3 [37] for the 
current and two future climate periods (2050 and 2070) 
(Figure 5). In addition to providing the details of invasion 
at the district level, we have also analyzed the elevation 
migration of parthenium weed to better determine the 
impact of climate change in this mountainous country. 
Our study will help direct future management plans and 
inform policy makers on the need for early detection and 
management in newly invaded areas. 


Area (Km’) 
PNW 
o888 


RCP 8.5 2050 
Loss Stable Gain 


260 Kilometers 


t t t 1 


Figure 5. Change in the distribution of P. hysterophorus between current and future climate scenarios. The maps were pre- 
pared in ArcGIS 10.4. Note: (i) Gain, (@®) loss, (@) stable, (~~) Not suitable, (>) Districts. 
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Besides climatic variables, there are other factors such as 
dispersal, demography and species interactions that de- 
termine species distribution [77]. However, information 
on such factors is not easily available and including them 
increases the complexity of the modelling process [78]. 
Additionally, in mountainous country like Bhutan where 
elevation ranges from 150 m asl to more than 7,000 m asl 
with steep slopes, topography can play a very critical role 
in creating barriers to the dispersal of species. Therefore, 
we included slope as one of the predictor variables in the 
modelling process. However, we did not find any influ- 
ence of slope in parthenium weed distribution in Bhutan, 
which did not agree with the findings of [23] in Nepal, 
where slope was found to have contributed the most. This 
was attributed to the choice of model algorithm where we 
have used an ensemble modelling approach instead of the 
single model modelling approach. If individual model was 
considered, MaxEnt did give slope as the most influencing 
variable (Table 2). 


Climate change has expanded the ranges of species 
[31,79]. Like other studies in the region [31,44,80,81], 
we have also demonstrated the impact of climate change 
and showed that the land suitability for parthenium weed 
invasion significantly increased under future climatic 
conditions. The percentage increase ranged from 7.39 to 
13.95% of Bhutan’s total area, or a 161 to 393% increase 
of the current predicted potential suitability. These find- 
ings are in sharp contrast to those of [12] where they 
have reported suitable land for parthenium weed inva- 
sion to contract in 2070 by about 4.15% as compared to 
the present 8.36%. Again, this may be attributed to the 
modelling method and choice of algorithm used, as they 
only used the MaxEnt program. Under current climatic 
conditions, out of 20 districts, 17 districts showed suit- 
ability to parthenium weed invasion. In the future, be- 
cause of climate change, two additional districts (Thim- 
phu and Paro) become suitable. Bumthang remained free 
from invasion under all the climatic conditions tested. 
These findings were also validated with the present day 
occurrence points and with the observations of agricul- 
ture extension officials working in those districts. Occur- 
rence points were collected only from the western and 
southern districts of Bhutan but the prediction showed 
invasion suitability in eastern and central districts where 
occurrence data were not collected such as Lhuentse, 
Mongar, Trashigang, Trongsa and Zhemgang. Additional- 
ly, [28] reported parthenium weed in Mongar, Trashigang 
and Trongsa while the first author detected five plants of 
parthenium weed in Thimphu around his residence area 
(27.51396°N, 89.6404°E, 2, 400 m asl) (model predicted 
suitability only in RCP8.5 2050 and 2070). 

Among the nine predictor variables used in the modelling, 
the most important ones that explained the species’ en- 
vironmental requirements best were the mean tempera- 


ture of the wettest quarter (Bhutan_Bio8), precipitation 
seasonality (Bhutan_Bio15), precipitation within the 
warmest quarter (Bhutan_Bio18) and the min tempera- 
ture of the coldest month (Bhutan_Bio6). The suitabili- 
ty of parthenium weed increased with increasing mean 
temperature of the wettest quarter (Bhutan_Bio8) and 
with precipitation seasonality (Bhutan_Bio15) while it 
decreased with increasing precipitation of the warmest 
quarter (Bhutan_Bio18) and the min temperature of the 
coldest month (Bhutan_Bio6). The low temperature of the 
mountains is the main limiting factor preventing species 
dispersal from lowlands into the mountains [82,83]. Ac- 
cording to the observations of [84] and [85] in Australia, 
the distribution of parthenium weed was limited to the 
temperatures between 5°C and 40°C. One of the major de- 
terminants of the parthenium weed distribution is mois- 
ture availability [86,87]. However, climate change will 
create climatically suitable niches in higher elevations 
through increased temperature and parthenium weed 
will thrive as a result, and as predicted by our model. Cli- 
mate projections by the National Center for Hydrology 
and Meteorology (NCHM) in Bhutan have predicted an 
increased surface temperature of between 0.8 to 1.6°C 
under RCP4.5 and between 0.8 to >3.2°C under RCP 8.5 
during the period 2021 to 2050 [42]. We have predicted 
a maximum upward elevation shift to 2,931 m asl; about 
753 m asl increase under RCP8.5 and by 2070. Presence 
of parthenium weed at 2,600 m asl was reported in Bhu- 
tan [13], while in Ethiopia, the highest altitude record of 
infestation was at 2,627 m [88]. The upslope movement 
of parthenium weed would mean highland pastures com- 
ing under threat from invasion and making such pastures 
unavailable for livestock production, especially for high 
altitude dwelling animals like yaks. 


Besides climate change, anthropogenic disturbance also 
drive the plant invasion in mountain ecosystems [89,83]. 
Bhutan has undergone unprecedented socio-economic 
development since its transition to democracy in 2008. 
Road network expansion is one of the main pledges of 
the political parties. Even those areas where construction 
was once restricted (e.g., Sakteng Wildlife Sanctuary) is 
now connected by roads. As roads are considered to be 
the major conduits through which parthenium weed is 
dispersed [90,88], roads will invariably allow partheni- 
um weed to invade not only in lowland regions but also 
at higher elevations in Bhutan. Anthropogenic distur- 
bance creating suitable habitat for parthenium weed may 
already be occurring as seen from the location of occur- 
rence points collected. The highest infestation and occur- 
rences were collected along the roads where one of the 
biggest hydropower projects in the country, Punatshang- 
chu Hydro-Power Project II is currently being built. It may 
be conjectured here that excavations and extensive move- 
ment of vehicles carrying hydro-electricity equipment 
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might have facilitated the spread of parthenium weed. 
It is more challenging to control invasive weeds in the 
mountainous areas than in adjoining lowland areas [91]. 
Therefore, early detection and the prevention of parthe- 
nium weed population establishment will be very critical. 


Bhutan is renowned for its very strong environmental 
protection policy, high density of plants and animals and 
the constitution necessitates that the country maintain 
60% of the land under forest cover in perpetuity [92]. 
Given the unfettered support for environmental protec- 
tion, the country’s mountainous terrain with steep slopes, 
high elevations, and dense forest may all contribute to 
preventing encroachment of invasive species. Howev- 
er, as stated, anthropogenic disturbances are increasing 
as political parties try to fulfill their pledges, and as the 
country strives to alleviate poverty and secure food for all. 
Further, Bhutan’s dependence on agriculture for econo- 
my and livelihood while only about 7% of the country’s 
area is arable, biological invasion may have a significant 
impact on food security and food production. Farmers 
have already reported impact of biological invasions in 
Bhutan [20]. Thirty horses have succumbed to Crofton 
weed consumption [19]. In Royal Manas National Park, 
61 ha of forest that was cleared to convert to grassland for 
wildlife has been invaded by Crofton weed [14]. Further- 
more, Bhutan is highly vulnerable to climate change [42]. 
While in-depth studies on the impacts of IAPS in Bhutan 
has not been undertaken, inventories are being initiated 
[10,93,16]. In this present study we have modelled the 
distribution of one of the most invasive weeds, partheni- 
um weed in Bhutan under the current and future climates 
and have shown that it will become highly invasive as a 
result of climate change and with the species even moving 
into higher elevation zones. This study should encourage 
land managers into undertaking early detection, mon- 
itoring and developing management plans for parthe- 
nium weed in Bhutan. One of the management options 
that could be explored is the use of biological agents, an 
approach that is in line with Bhutan’s policy of restricted 
chemical use for environmental protection. The presences 
of two biological agents, the Mexican beetle (Zygogramma 
bicolorata Pallister, 1953) and the winter rust (Puccinia 
abrupta var partheniicola (H.SJacks.) Parmelee) have 
already arrived in Bhutan [94]. These agents have been 
introduced and have successfully established in Australia 
[95] and India [96] and are an important component of 
those countries management plans of parthenium weed. 


Conclusion 


We have modelled the impact of climate change on the 
distribution of one of the most invasive weeds, partheni- 
um weed in Bhutan. Through our study, we have shown 
that under the future climate period (2050 and 2070) the 
suitable habitat areas for parthenium weed will expand 
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by as much as five times the current distribution, and 
come to infest all districts of Bhutan, except for the high 
mountainous district of Bumthang. Under the future cli- 
mate, the highest land suitability for parthenium weed is 
predicted under RCP4.5 2050 period with about 5,419.69 
km’. Generally, the districts in west and south showed 
more suitability than the districts in the east and central. 
The highest elevation of suitability was predicted at 2,931 
m asl; an upward shift of about 753 m by 2070. These 
observations suggest that many of the country’s native 
species will be at risk from displacement by partheni- 
um weed and might even be pushed to extinction. Food 
security and food production in Bhutan will also be con- 
strained as parthenium weed commonly invades agricul- 
tural areas. Therefore, early detection and eradication 
plan of new parthenium weed populations will be very 
important as currently, only about 2.83% of the country’s 
area is predicted to have suitability to parthenium weed 
invasion, which appears to be quite manageable. Early 
detection and rapid response should remain our priority 
since delay in control programs will be very costly with 
eradication impossible in future. 
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